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Q_|. Both analytic and numerical radiative transfer models applied to high spec- 

Q . tral resolution CS and N2H"^ data give insight into the evolutionary state of 

il3 . L1551 MC. This recently discovered pre-protostellar core in L1551 appears to 

d • be in the early stages of dynamical evolution. Line-of-sight infall velocities of 

. > O.lkms"^ are needed in the outer regions of L1551 MC to adequately fit the 

^ . data. This translates to an accretion rate of ~ 10 Mq yr , uncertain to within 

a factor of 5 owing to unknown geometry. The observed dynamics are not due 
to spherically symmetric gravitational collapse and are not consistent with the 
standard model of low-mass star formation. The widespread, fairly uniform CS 
line asymmetries are more consistent with planar infall. There is modest evi- 
dence for chemical depletion in the radial profiles of CS and C^^O suggesting 
that L1551 MC is also chemically young. The models are not very sensitive to 
chemical evolution. L1551 MC lies within a quiescent region of L1551 and is 
evidence for continued star formation in this evolved cloud. 

Subject headings: ISM: clouds — ISM: individual(L1551) — stars: formation — 
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1. Introduction 

Gravitationally bound molecular cores with no embedded infrared sources, referred to as 
pre-protostellar cores (PPCs), offer the unique opportunity to study the conditions preceding 
stellar birth. To date there have only been a few PPCs studied in depth leaving many 
unanswered questions regarding the kinematics and chemistry during this epoch of star 
formation. 

Inward gas motions are necessary for the formation of stars and the kinematic sig- 
natures of these motions are expected to be present in molecular line data. Blue-peaked 
asymmetric line profiles commonly seen in optically thick tracers of molecular cores can be 
accurately modeled as self-absorbed infalling gas (Evans 1999; Myers et al. 1996). A high 
ratio of blue-peaked to red-peaked profiles in selected cores is further evidence that these 
spectral signatures are tracing global infall in some PPCs (Lee et al. 2001). These "infall 
asymmetries" imply infall velocities of up to ~ 0.1 kms^^ extending out to ~ 0.1 pc (Tafalla 
et al. 1998; Lee et al. 2001) inconsistent with the ambipolar diffusion paradigm of isolated 
low-mass star formation (Shu, Adams & Lizano 1987; Ciolek & Mouschovias 1995). 

The observed decrement in abundance of certain molecules such as CO and CS toward 
the centers of cores (Kuiper et al. 1996; Willacy et al. 1998; Caselli et al. 1999; Tafalla et al. 
2002) can be explained by the expected depletion of these species onto dust grains (Watson 
& Salpeter 1972; Bergin & Langer 1997). In some PPCs the abundance profile is best fit 
using a model with a central region of complete depletion (e.g., Tafalla et al. 2004). Models 
have shown that depletion becomes more severe as a core collapses suggesting that chemical 
depletion may be useful in determining the age of PPCs (Rawlings et al. 1992; Bergin & 
Langer 1997; Lee et al. 2004; Aikawa et al. 2005). However, detailed molecular fine studies 
have shown that chemistry is not reliable as a sole evolutionary indicator (Lee et al. 2003) 
and interpretations should be made on a case by case basis (Rawlings & Yates 2001). 

Analytic models of asymmetric molecular line profiles (e.g. Myers et al. 1996) provide 
quick solutions by imposing functional forms for molecular excitation and chemistry along 
the line of sight. These models have been shown to reliably characterize simulated data 
despite using the simplest available assumptions (De Vries & Myers 2005). Numerical models 
typically contain less assumptions and are hence more fiexible, but are computationally much 
more costly. 

L1551 MC is a PPC that sits in a relatively isolated region to the northwest of IRS5 in 
the well-known L1551 dark cloud (Swift et al. 2005, hereafter Paper I). It is expected to be 
gravitationally unstable and hints of infall have been seen in position switched CS spectra. 
The differentiation seen between the NH3 and C2S emission suggests chemical evolution. 
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Past studies of the L1551 cloud (e.g., Snell 1981; Devine, Reipurth & Bally 1999) provide 
a well-defined context for L1551 MC. Following the discovery and primary characterization 
of this starless molecular core (Paper I), this article investigates dynamical and chemical 
signatures from new molecular line data in an attempt to understand its state of evolution 
toward stellar birth. 

In § 2 we outline the molecular line observations used with the observations of Paper I 
to analyze L1551 MC. The models used to interpret these data are described in §3 with 
the results presented in § 4. Further consideration of these results in § 5 reveal the main 
conclusions of this paper which are then summarized in § 6. 

2. Observations and Reductions 

We used the Five College Radio Astronomy Observatory's 14 m telescope on the nights 
of 2005 January 7 and 8, and for partial nights on February 8, 9 and 24 to map L1551 MC 
in N2H^and CS simultaneously. The weather was variable on the nights of January 7 and 
8 with system temperatures varying from ~ 220 K to ~ 375 K. The partial nights were 
substantially better showing system temperatures from 130 K to 150 K. 

The observing frequencies in the first and second IF of the dual channel correlator (DCC) 
were set to 93.176258 GHz and 97.980953 GHz, corresponding to the isolated component of 
the N2H"''(1 — 0) and the CS(2 — 1) rotational transitions (Lee et al. 2001) respectively, for 
our maps. The beam at these observing frequencies has a FWHM of ~ 45" and a main 
beam efficiency rj^ ~ 0.5. The highest spectral resolution mode of the DCC was used giving 
a 25 MHz band with 1024 channels resulting in ~ 0.08 km s~^ velocity resolution. 

A 6'x6' patch of sky was mapped centered on L1551 MC, 4^^31"09".9, +18°12'41" (J2000.0), 
using the on-the-fly observing technique. Our final maps are a composite of ~ 20 individ- 
ual maps created by scanning alternatively in the right ascension and declination directions 
combined by weighted mean. These maps are highly oversampled given the 12" row spacing, 
60" s""^ scan speed and the rotation of the SEQUOIA 32 element receiver array with par- 
alactic angle. The data are calibrated with an off scan every two rows and a total power 
calibration scan every four rows. These data were reduced, griddcd and output to FITS files 
using the OTFTOOL utility available at the observatory. The RMS noise level in our final CS 
and N2H"'" maps is ~ 26 mK. All further analysis was performed using the IDL software of 
Research Systems Inc. 
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3. Models 

We use a non-linear least squares minimization based on the Levenberg-Marquardt 
technique to find the best fits to our spectral line data according to two distinctly different 
radiative transfer models. The analytic model fits the asymmetric CS line profiles while the 
family of numerical models use information from the optically thin N2H+ emission to fit both 
the CS and the CO isotopologue data. 

3.1. Analytical Model 

The HILL5 model presented in De Vries & Myers (2005) assumes a symmetric 1-dimensional 
structure for all parameters. The central excitation temperature, Tc, tapers linearly to the 
cloud edge where it has a value of T^, — 2.73 K. Gas approaches the center of this model 
system at a constant velocity, relative to the systemic velocity, vlsr- Including the 
central optical depth, tc, and the velocity dispersion of the gas, o", there are a total of 5 free 
parameters. The hill5 model is the preferred model in De Vries & Myers (2005) giving 
estimates of the infall velocity in simulated observations good to 0.02 km s~^. 

3.2. Numerical Models 

The "rt" program introduced in Dickel & Auer (1994) (also see Williams et al. 2004) 
computes the full, non-equilibrium radiative transfer for a spherically symmetric mass. This 
FORTRAN code takes the density, kinematics, kinetic temperature and molecular abundance 
ratio as a function of core radius and solves for the level populations using a lambda-iteration 
technique. 

Table 1 outlines the fixed and free parameters which define our family of numerical 
models and Table 2 summarizes the individual models described in the following subsections. 
The density of molecular hydrogen is represented by n. The density of a molecular species 
is denoted by n„toi and is related to the total density by the conversion factor, X„ioi- The 
infall velocity, is equal to vq at the minimum radius of the model and follows a power 
law form with index p. The model core is isothermal with a constant turbulent velocity 
width /\.Vturh- Quantities are evaluated at 40 logarithmically spaced radii from r^j„ to r^aa;- 
The fixed parameters, except for the radii, are derived from the NH3 data of Paper I. The 
maximum radius represents the extent of the CS emission. The minimum radius was chosen 
in consideration of a timely convergence and our results are not sensitive to its exact value. 
All models have as free parameters no, vq, and X^oi- 
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Each numerical model comprises a substantial parameter space to explore. The density 
and infall velocities are therefore constrained to 10^ < no < lO^cm"^ and < Wj„ < 
0.5 km s~^ in accord with the results of Paper I. For each model the best fit was determined by 
repeating the least-squares minimization procedure beginning at randomly chosen parameter 
values until the global minimum in -y^ space was converged upon. 

These models were first used to fit the CS(2 — 1) line profiles and then modified to fit 
for the conversion factors for ^^CO and C^^O given the best fit parameters from the CS fits. 
The kinetic temperature and turbulent velocity width are also free parameters in these fits 
to account for the possibility that the CO isotopologucs arc tracing different environments 
within L1551 MC. Upper limits of 12 K (Snell 1981) and 0.5 km s^^ were enforced for Tk and 
^Vturh respectively. 

The systemic velocity for all the numerical models is obtained from fits to the full 
N2H"'"(1 — 0) transition with the hyperfine fitting routine developed in Paper 1. This routine 
returns the optical depth, velocity centroid, velocity width, and excitation temperature for 
a given spectrum. 

The normalized root mean square of the error is used as our goodness of fit indicator 
defined by 



where a is the RMS noise level in the spectrum, iV is the number of data points used in the 
fit, Di are the data points and yi are the model points. Values of RMSE ~ 1 are expected 
for successful fits. 

There are six total numerical models denoted by an identifier and a Roman numeral 
(Table 2). Each pair of models are identical except for the abundance profile which is either 
constant throughout or takes the form 




(1) 
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which is referred to as a "drop" abundance profile. The models are further grouped by the 
assumed velocity structure into two main categories, accretion and collapse. 
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3.2.1. "Accretion" Models 

The first four models are referred to as "accretion" models because they allow for the 
possibility of high velocity gas in the outer regions of the core. Model ACV-I is like the 
HILL5 analytic model, but the radiative transfer is solved for numerically. Therefore, the 
infall velocity from ACV-l can be compared to the results of the hill5 model for consistency. 
Model ACV-ll is used to test whether or not having a depletion zone significantly alters the 
derived infall velocity. 

Models A-l and ll have a velocity profile that increases as the square root of the core 
radius and represent scenarios in which gas is accreting onto the core from large radii. These 
models, in comparison with the others, test if the highest velocity infalling gas is in the 
exterior regions of the core. 

3.2.2. "Collapse" Models 

Models GC-l and ll have velocity structures chosen to represent inside-out collapse, as 
in Shu (1977). Given the complexity of L1551 MC from observations a detailed comparison 
with theoretical models would be difficult. However, free-fall collapse is generally expected to 
produce a decreasing velocity profile with radius from the center of the core to the rarefaction 
front (e.g., Fatuzzo et al. 2004; Myers 2005). Therefore these models test if the motions 
implied by the CS emission are consistent with the motions expected from gravitational 
collapse directly associated with the formation of a protostar. 

4. Results 

Figure 1 shows the distribution of velocity integrated N2H+ in greyscale overlaid with 
independent CS(2 — 1) spectra at their respective plane-of-sky locations. The vertical dotted 
lines shown in selected spectra are the systemic velocities from the hyperfine fits to the 
N2H+ profiles. Most CS line profiles are asymmetric having a peak blueward of the systemic 
velocity and a smaller red peak or shoulder. 

4.1. Statistics of the CS Line Asymmetry 

The normalized velocity difference introduced in Mardones et al. (1997) is a useful 
quantity to characterize spectral asymmetries. Using CS and N2H+ as optically thick and 
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thin molecular species, we define 

^^^^ = ~7^/ 

where Vcs is the velocity centroid of the brightest peak in the CS hne profile and the velocity 
centroid and intrinsic FWHM of the N2H+ emission, Vn2H+ and ^Vn^h+, are obtained from 
hyperfine fitting. 

The N2H+ emission is less extended than the CS emission, so it is the 31 beam-sampled 
spectra with well-determined centroids ((Jv^^^ < 0.1 Ay;v2ij+) that define our sample. Figure 2 
shows a histogram of the 5Vcs values. The cross-hatched area corresponds to 5Vcs values 
that are greater than 5 times the estimated error [as) obtained from standard propagation 
of error. The values for as range from 0.03 to 0.4 with a median of 0.10. The mean value 
{5Vcs) = -0.50. 

The blue excess is defined as £^ = (iV_ — N+)/N (Lee et al. 2001), where iV_ is the 
number of positions with 5Vcs ^ — Scr^, iV+ is the number of positions with 5Vcs ^ 5o"5 and 
here = 31. The blue excess of L1551 MC is high, E = 0.39, meaning a significant fraction 
of CS lines have a blue peak. 

The P value of the student t-test measures of the probability that a given sample of 
values was drawn from a normal distribution with a mean of zero. For the 31 5Vcs values 
in L1551 MC, P is less than 0.5%. Together with a high E value, this means that the 
asymmetric CS profiles are not random but are strongly biased toward having a blue peak. 
In relation to the cores studied in Lee et al. (2001), L1551 MC would rank as a bona fide 
infall candidate. 



4.2. Model Fits 

The 18 spectra with I^Vcsl > 5o"5 are fit according to the HILL5 model. Three of these 
18 spectra return positive (outflowing) gas velocities and correspond precisely to the three 
spectra that show positive 5Vcs values in Figure 2. Of the 15 spectra that show negative 
(infalling) gas velocities, the mean = 0.13kms~^, roughly two thirds the sound speed 
in the medium. 

Averages of our molecular line data within the central 100" of L1551 MC (dashed circle 
in Figure 1) yield high signal-to-noise composite spectra in CS, N2H"'", ^^CO, and C-'^^O. The 
best fit to the CS composite spectrum using the hill5 model is outlined in Table 3 and is 
presented as the dotted line on the CS spectrum of Figure 3. The fit is good overall, but 
slightly overestimates the emission in the red wing of the profile. 
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The results from the numerical fits are summarized in Table 4. The best fit to the CS 
composite spectrum was achieved with model A-ll with an RMSE value of 0.77 and is shown 
as the dashed line on the CS spectrum of Figure 3. The constant velocity models, ACV-l 
and II fit the data fairly well and are in agreement with the HILL5 results, while the collapse 
models, GC-l and ll, produce the worst agreement with the data. 

5. Discussion 
5.1. Dynamics 

There is a surplus of blue peaked asymmetric CS line profiles in L1551 MC successfully 
modeled as self- absorption from globally infalling gas (§4). In both our analytic and numer- 
ical models, gas velocities > O.lkms"^ are needed to account for the observed CS profiles. 
These motions are not associated with any of the known outflows in the L1551 cloud, in- 
cluding those from IRS5 (Snell, Loren & Plambeck 1980; Moriarty-Schieven & Snell 1988; 
Moriarty-Schieven & Wannier 1991), L1551 NE (Pound & Bally 1991; Moriarty-Schieven et 
al. 1995; Devine, Reipurth & Bally 1999) and XZ/HL Tau (Welch et al. 2000). A closer look 
at our results reveals insights into the nature of these observed dynamics. 

5.1.1. "Collapse" vs. "Accretion" 

Figure 4 shows a direct comparison between the best "collapse" model flt, GC-l, and the 
best "accretion" model fit, A-ll. The GC models produce the least consistent fits to our data. 
In these models the absorbing gas is in the outer regions of the core where the excitation 
temperatures are lowest and where gas is moving at relatively slow velocities. Therefore, 
these models produce centralized self-absorption with a double peaked profile instead of the 
observed red shoulder. Higher values of vq push the absorption feature slightly further to 
the red and this is why the best fits for the gravitational collapse models are pegged against 
the Vq — 0.5 km s"^ limit set by our constraint on highly supersonic gas (§3.2). 

The constant velocity accretion models, ACV-I and II, do substantially better than the 
inside-out collapse models allowing for absorption and emission of high velocity gas. But 
the monotonically increasing infall velocity of the A-l and ll accretion models works together 
with the decreasing CS excitation temperature as a function of radius to carve a smooth red 
shoulder in the line profile most consistent with our data. A marginally better fit to the data 
is achieved with the drop abundance profile of A-ll and we use this as the preferred model. 
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5.1.2. High Velocity Gas at Large Core Radii 

The observed self-absorption in the CS hne profiles forces our numerical model fits into 
the optically thick regime. It is therefore the dynamics, rather than the chemistry or mass 
distribution, that most significantly affects the emergent line profiles. Our fits show that 
infalling gas is needed in the outer regions of the core to reproduce the observed CS line 
profiles. The hill5 model fits to the ensemble of CS spectra also show high infall velocities 
at large projected core radii (see §5.1.3). 

According to the standard model of low-mass star formation (Shu, Adams & Lizano 
1987), a core collapses via ambipolar diffusion until it becomes supercritical and then un- 
dergoes an inside-out collapse. The high infall velocities at large core radii indicated by 
our models are comparable to other closely studied cores where the dynamics were found 
to be inconsistent with the standard model (e.g. Tafalla et al. 1998; Lee et al. 2001). Prom 
our models it is clear that the gas motions observed across L1551 MC do not represent in- 
falling gas due to gravitational collapse and have magnitudes too large to be associated with 
ambipolar diffusion. 

The high infall velocities observed over a large angular extent around L1551 MC are 
more consistent with planar accretion. Planar accretion would produce this observational 
signature for a wide range of inclination angles. The magnetic field around L1551 has been 
measured to be fairly uniform (Vrba, Strom, & Strom 1976) and a roughly constant field 
could facilitate a planar accretion scenario. However, measurements of the magnetic field 
geometry in L1551 MC are essential to clarify this situation. 

Prom the analytic and numerical modeling of N2H"'" and CS spectra in L1551 MC it 
appears that gas is falling toward the mid-plane of the core from large radii. The observed 
dynamics are not directly related to the formation of a star, but do not preclude gravitational 
collapse scenarios (e.g., Patuzzo et al. 2004). 

5.1.3. Mass Accretion 

High infall velocities are derived from the HILL5 model for a vast majority of the CS 
spectra which span ~6', or 0.24 pc at the distance of Taurus. This loosely defines a radius 
at which material is infalling, Rm ~ 0.1 pc. 

The mass infall rate is estimated as the amount of material passing through a spherical 
surface with radius Rin, = AjiRfi^pVin. Here we use a density of 1.2 x lO^cm"^ and 
an infall velocity of 0.02 km s^^ derived from the model A-il best fit parameters and the 
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density and velocity relations from Table 1 evaluated at Rin- The derived infall rate Mj^ ~ 
10~^ Mq yr~^. Values of the average density (Paper I) and average infall velocity (§ 4.2) yield 
an accretion rate of the same order of magnitude as does planar infall at a constant velocity 
of O.lkms^^. Due to the unknown geometry and symmetry in this problem the accretion 
rate cannot be constrained within a factor of ~ 5. 



5.2. Chemical Structure 

5.2.1. Comparisons of Molecular Species 

A fit to the mean spectrum within a half power contour of the N2H^ velocity integrated 
emission gives a central velocity of 6.71 km s~^, an observed line width of 0.28 km s~^ trans- 
lating into a non-thermal width of 0.25 km s~^ for a kinetic temperature of 9K and a total 
optical depth of 4.03 (r = 0.4 for the isolated component). This is very similar to the NH3 
emission with vlsr — 6.72 km s~^ and Avnt — 0.2kms~^ (Paper I). The sky distribution of 
N2H+ also follows the NH3 closely with a major and minor FWHM of 2.'9 and 1.'3 oriented 
at a position angle of 139° compared to the NH3 values of 2.'3, l.'l and 133°. The emis- 
sion from both molecular species peak at the same sky coordinates and coincide precisely 
with a peak in the sub- mm dust continuum (Schieven 2004). The N2H+ velocity gradient 
across the core is 1.3kms~^ pc~^ at a position angle of 227° compared to the NH3 values of 
1.2kms~^pc~^ and 224°. Therefore it seems the N2H''" and NH3 emission trace the same 
high-density interior gas in L1551 MC. 

Figure 5 compares the radial profiles of NH3, N2H+, CS, and C^^O. The NH3 and N2H"'" 
profiles are fit using Equation 1 and the CS and C^^O profiles are fit using Equation 2 of 
Paper I. The radial distribution of N2H+ is slightly flatter than the NH3 which may be 
due to a slight overabundance of NH3 in the innermost regions as seen in other Taurus 
cores (e.g. Tafalla et al. 2004). The CS and C^^O profiles are significantly flatter than the 
NH3 and N2H"'" suggesting that the relative abundance of these species is higher at larger 
radii. While this may be due in part to the optical depth of CS, the C^^O is optically thin in 
this region (see §5.2.2). Figure 6 shows that the CS has a much broader distribution on the 
sky compared to the N2H"'" despite the similarity in critical density values for the transitions. 
A depression in the CS emission ~ 1' northwest of the peak N2H"'" and NH3 emission is also 
seen in Figure 6 that coincides with a deflcit in C^^O emission. This shght depression is 
bordered on the southwestern edge by the C2S emission and may signify a depletion zone. 

It is not clear whether or not significant depletion has occurred in the inner regions 
of L1551 MC since chemical effects cannot be sufficiently isolated with our models. How- 



- 11 - 



ever, the relative deficit of emission from both C^^O and CS at small radii indicate some 
depletion of these species has occurred presumably due to chemical evolution related to core 
contraction (e.g. Lee et al. 2004). 

5.2.2. CO Isotopologue Emission 

Figure 3 shows the best fit to the composite CO isotopologue spectra and Table 5 
outlines the fit results. The best fit turbulent velocities for our CO isotopologues are a 
factor of 1.5 higher than the CS and N2H"'" fits. These larger values are likely due to the 
CO isotopologues tracing a larger volume along the line of sight. The CO isotopologues are 
tracing cold gas with the ^^CO and C^^O being optically thick and thin respectively. 

The abundance of both isotopologues derived from our models are consistent with ob- 
servations (e.g., Frerking et al. 1982: Lada et al. 1994) given the stated uncertainties. The 
best fits for Xwqq and Xq^sq arc 1.1 x 10^^ and 1.1 x 10^'' respectively. For all our models 
the ratio of these abundances most consistent with the observed CO isotopologue line inten- 
sities is ~ 10; a factor of 2 higher than what would be estimated from terrestrial abundance 
ratios (Pickett ct al. 1998). 

The centroid of the C^^O line agrees well with the centroid of the N2H^ (see Figure 3) 
in accord with the results from Lee et al. (2001). The shape of the ^^CO line profile, 
slightly broadened by optical depth, is fit well with model A-ll, but the central velocity is 
clearly shifted. This shift accounts for the anomalously high RMSE values in Table 5. A 
cross correlation between the model A-ll fit and the observed spectrum peaks at a lag of 
7.1kHz corresponding to a velocity offset of 0.02 km s~^ at this transition frequency. This 
is not accounted for in the stated measurement errors for the published frequency equal to 
5.1kHz (Pickett ct al. 1998). If this offset were due to an error in the chosen observing 
frequency, a best estimate of the true frequency would be 110.2013470 GHz. 

5.3. Is the C2S Emission Tracing the CS Self-Absorbing Gas? 

It is likely that CS and C2S share the same chemical pathway in the interstellar medium (e.g., 
Suzuki et al. 1992). Yet the critical density of CS(2 — 1) is an order of magnitude higher 

than C2S(32 — > 2i) (Wolkovitch ct al. 1997). It is therefore possible for these species to 
exist in the same environment in different excitation states. The C2S emission in L1551 MC 
is shifted to the red of the systemic velocity (see Paper I) while the CS spectrum shows 
redshifted self-absorption. 
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Figure 7 shows composite spectra of the CS and C2S hnes averaged within a half maxi- 
mum power contour of the C2S emission (see Figure 6). The composite CS spectrum shows 
a blue peaked, asymmetric profile. The best fit to this spectrum using the hill5 model 
returns an infall velocity of 0.14kms~^, nearly identical to the velocity centroid of the C2S 
emission. 

This supports the idea that the C2S emission is tracing the CS absorption layer, but it 
still remains unclear why there is no blue component to the C2S profile. Chemistry is likely 
responsible. 

6. Conclusions 

High spectral resolution maps of CS and N2H+ taken with the FCRAO 14 m telescope are 
fit with two distinctly different radiative transfer models to gain insight into the dynamical 
and chemical evolution of the recently discovered PPC in Taurus, L1551 MC. 

The analytic and numerical models both suggest that L1551 MC is in the early stages 
of dynamical evolution. Infall velocities of > O.lkms"^ are needed in the outer regions of 
L1551 MC to fit our data adequately. This translates to an accretion rate of ~ 10~^ Mq yr"-*-, 
uncertain to within a factor of 5. The observed dynamics are not due to gravitational collapse 
and are not consistent with the standard picture for low mass star formation. The extent 
and shape of the asymmetric CS line profiles are more consistent with planar accretion. 

There is modest evidence for chemical depletion near the highest density regions of 
L1551 MC in CS and C^^O. This suggests that L1551 MC is chemically young. Our models 
are not very sensitive to chemical evolution and do not shed light on this issue. 

Fits to our CO isotopologue data yield a ^'^CO to C^^O relative abundance ratio of ~ 10. 
These fits also reveal the ^'^CO line centroid to be shifted relative to the other molecular 
species by 7.1 kHz or 0.02 km s^^. If this is due to an incorrect observing frequency, the best 
estimate for the true transition frequency is 110.2013470 GHz. 

The C2S emission seen in Paper I may be tracing the infalling CS gas that is producing 
the observed self- absorption. This would explain the apparent velocity offset of C2S with 
respect to the optically thin, high density tracers. However, it is unclear why there is no 
blue component to the C2S line profile. 

The gas dynamics and distributions of molecular species in L1551 MC are distinct from 
the complex processes related to the known young stellar objects within the greater L1551 
region and are evidence for continued star formation in this evolved cloud. 
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Fig. 1. — Velocity integrated N2H+ emission from L1551 MC in greyscale with independent 
spectra of CS(2 — 1) overlaid at their respective positions. The dotted contour represents the 
50% peak N2H+ emission contour. The dashed circle denotes the region in which composite 
spectra were constructed for N2H+, CS, ^^CO, and C^^O. The vertical dotted lines in selected 
spectra represent the systemic velocity along the line-of-sight determined from fits to the 
N2H+ spectra at the respective locations. 
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Fig. 2. — Histogram of SVcs values for all spectra with reliable velocity centroids determined 
from fits to the N2H+ data. The shaded region corresponds to the histogram of spectra with 
If^Vcsl > 5(75. See §4.1 for details. 
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Fig. 3. — N2H+, CS, C^^O, and ^^CO composite spectra in brightness temperature units. 
The scale of the C^^O and ^^CO spectra are reduced by a factor of 2 and 5 respectively. The 
model A-II fits to all spectra and hill5 fit to the CS spectrum are shown as dashed and 
dotted lines respectively. The thin vertical line signifies the velocity centroid of the N2H+ 
line. 
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Fig. 4. — Comparison between the best fit "collapse" model, GC-l, to the best fit "accretion" 
model, A-ll. The solid line denotes the data while the model fits are represented with the 
broken line. See §5.1.1 for details. 
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Fig. 5. — Profiles of four different molecular species averaged in elliptical annuli around 
L1551 MC. For each species the flux is normalized to the peak value in the proflle. Equation 
1 (for NH3 and N2H+) and Equation 2 (for CS and C^^O) from Paper I are used to flt the 
proflles. 
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Fig. 6. — Velocity integrated CS in greyscale overlaid with contours of N2H+ (black contours) 
and C2S (white contours; Paper I) at 30, 50 70 and 90% peak flux. The beam size for the 
CS and N2H+ data is shown at the bottom right. The black triangle denotes the position of 
HH265. 
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Fig. 7. — CS and C2S spectra averaged over the half maximum contour of C2S emission in 
L1551 MC. The vertical dotted line shows the central velocity of the infalling layer determined 
from the hill5 fit to the CS line profile and is seen to be nearly identical to the velocity 
centroid of the C2S line profile fit by a single Gaussian component. 
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Table 1. Overview of Numerical Cloud Models 



Model Parameter Value 

n(r) no/ [1 + ir/r^"] cm~^ 

rimoi {r) X^oi {r^ n{r) cm'^ 

v{r) vo {r/rrainY kms"^ 

/!^Vturb 0.2 km s"-*^ 

Tk 9K 

Tmin 67 AU 

Tmax 0.2 pC 

Tc 0.05 pc 

a 2.2 



'^Xjnoii''^) is either constant with radius or 
takes the form of Equation 2. 



Table 2. Family of Numerical Models 



Model Name Velocity Profile Abundance Profile Comments 



Accretion" Models 



ACV-I 


v{r) 


= constant 


X„ioi{r) 


= constant 


Numerical analogue to hill5 model. 


ACV-II 


v{r) 


— constant 




= "drop"^ 


Like model ACV-I but with central 












depletion. 


A-I 


v{r) 


oc r^/^ 


Xmoi{r) 


= constant 


High velocity gas at large core radii. 


A-II 


v{r) 


oc r^/^ 


Xmoi{r) 


= "drop"^ 


Like model A-l but with central 












depletion. 



"Collapse'^ Models 



GC-l v{r) oc r~^/^ Xmoi{r) — constant Inside-out collapse model with 

constant abundance fraction. 

GC-ll v{r) oc r~^/^ Xmoiif) — "drop"^ Inside-out collapse model with 

central depletion. 



^See Equation 2. 
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Table 3. hill5 Model Fit to Average CS Spectrum 

Model Parameter Best Fit Value 

rc 2.02 

(7 O.lSkms"^ 

Tc 4.57 K 

vlsr 6.72 km 

Vin 0.15 km 

RMSE*^ 2.03 



^The root mean square of the error 
as defined in Equation 1. 
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Table 4. Numerical Model Fits to CS Emission 



Name Uq Vq Xcs Tmax RMSE^ 

lO^cm-3 kms"^ 10-9 



■Accretion" Models 



ACV-I 


4.6 


0.16 


0.8 


2.8 


0.95 


ACV-II 


7.9 


0.15 


1.1 


5.1 


1.09 


A-I 


3.3 


0.015 


1.4 


3.0 


0.89 


A-II 


6.9 


0.001 


1.4 


5.2 


0.77 






"Collapse" 


Models 






GC-I 


4.2 


0.5 


1.1 


5.1 


2.38 


GC-II 


6.3 


0.5 


2.1 


8.9 


2.68 



^The root mean square of the error as defined in 
Equation 1. 
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Table 5. Numerical Model Fits to CO Isotopologue Emission 



Name 




Avt„r6[kms ^] 




^max 


RMSE^ 


ACV-I 


6.6(-7) 


0.39 


8.8 


3.1 


4.19 


ACV-II 


9.8(-7) 


0.39 


8.7 


4.9 


4.11 


A- 1 


8.6(-7) 


0.37 


9.0 


2.7 


4.47 


A-II 


l.l(-6) 


0.38 


8.8 


4.7 


4.18 


GC-I 


7.9(-7) 


0.42 


8.7 


3.6 


3.96 


GC-II 


1.3(-6) 


0.43 


8.7 


5.2 


3.96 




ACV-I 


6.9(-8) 


0.29 


9.0 


0.37 


0.96 


ACV-II 


1.0(-7) 


0.30 


9.2 


0.58 


0.95 


A-I 


9.2(-8) 


0.24 


10.3 


0.28 


1.02 


A-II 


l.l(-7) 


0.28 


12.0 


0.37 


0.99 


GC-I 


1.0(-7) 


0.34 


6.1 


0.97 


0.94 


GC-II 


1.7(-7) 


0.35 


6.1 


1.4 


0.94 



'^The root mean square of the error as defined in Equa- 
tion 1. 



